Leave-one-site out LICA

# Before ComBat
# all site (IC maps)
# Components were ordered by variance explaining.
all=readMat('./2023_11_14_lica_n140_ctareadti_B1000nn_rsfc.mat')

all_nocb = do.call(rbind,lapply(1:5, function(ss)all$Morig40[[7]][[ss]][[1]] %*% diag(as.vector(all$Morig40[[3]][[ss]][[1]]))))[,order(all$Morig40[[15]], decreasing=TRUE)]
colnames(all_nocb) = paste('noCombat.All.IC',stringr::str_pad(1:40,2, pad = "0"), sep='')
 

all_cb = do.call(rbind,lapply(1:5, function(ss)all$Morig40.combat[[7]][[ss]][[1]] %*% diag(as.vector(all$Morig40.combat[[3]][[ss]][[1]]))))[,order(all$Morig40.combat[[15]], decreasing=TRUE)]

colnames(all_cb) = paste('Combat,All.IC',stringr::str_pad(1:40,2, pad = "0"), sep='')

# leave one out site (IC maps)
# Components were ordered by variance explaining.

site=readMat('2023_11_14_lica_n140_ctareadti_B1000nn_rsfc_by_site.mat')
xmat_no_cb=do.call(cbind,lapply(1:5,function(j){
  tmpmat = do.call(rbind,lapply(1:5, function(ss)site$loso[[j]][[1]][[7]][[ss]][[1]] %*% diag(as.vector(site$loso[[j]][[1]][[3]][[ss]][[1]]))))[,order(site$loso[[j]][[1]][[15]], decreasing=TRUE)]
  colnames(tmpmat) = paste('Site',j,'.IC',stringr::str_pad(1:40,2, pad = "0"), sep='')
  return(tmpmat)}
  ))

xmat_cb=do.call(cbind,lapply(1:5,function(j){
  tmpmat = do.call(rbind,lapply(1:5, function(ss)site$loso.combat[[j]][[1]][[7]][[ss]][[1]] %*% diag(as.vector(site$loso.combat[[j]][[1]][[3]][[ss]][[1]]))))[,order(site$loso.combat[[j]][[1]][[15]], decreasing=TRUE)]
  colnames(tmpmat) = paste('Site',j,'.IC',stringr::str_pad(1:40,2, pad = "0"), sep='')
  return(tmpmat)}
  ))

# Find the best matching components to the all site
cor_no_cb = abs(cor(all_nocb, xmat_no_cb, method='pearson'))
cor_no_cb = cor_no_cb[, apply(cor_no_cb, 2, function(s)all(is.na(s)))==FALSE]
cor_no_cb = cor_no_cb[apply(cor_no_cb, 1, function(s)all(is.na(s)))==FALSE, ]
cor_no_cb2 = cor_no_cb[, apply(cor_no_cb, 2, function(s)all(abs(s<0.5)))==FALSE]
cor_no_cb2 = cor_no_cb2[apply(cor_no_cb2, 1, function(s)all(abs(s<0.5)))==FALSE,]

cor_cb = abs(cor(all_cb, xmat_cb, method='pearson'))
cor_cb = cor_cb[, apply(cor_cb, 2, function(s)all(is.na(s)))==FALSE]
cor_cb = cor_cb[apply(cor_cb, 1, function(s)all(is.na(s)))==FALSE, ]
cor_cb2 = cor_cb[, apply(cor_cb, 2, function(s)all(abs(s<0.5)))==FALSE]
cor_cb2 = cor_cb2[apply(cor_cb2, 1, function(s)all(abs(s<0.5)))==FALSE,]

#apply(cor_no_cb.tmp,2, which.max)
colist = c(rep("#FFFFE5",50),COL1('YlOrBr', 100))

drawcorplot <- function(mat, threshold=0.5){
  mat[abs(mat)<threshold]<-0
  mat = mat[apply(mat, 1, function(s)all(abs(s<threshold)))==FALSE,]
  mat = mat[,apply(mat, 2, function(s)all(abs(s<threshold)))==FALSE]
  
  corrplot(t(rerow(t(mat))),col = colist, col.lim=c(0,1),is.corr = FALSE)
}

Spatial pattern matching

no combat

drawcorplot(cor_no_cb2[,grep('Site1.', colnames(cor_no_cb2))], threshold=0.5)
## [1] "16 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."

drawcorplot(cor_no_cb2[,grep('Site2.', colnames(cor_no_cb2))])
## [1] "23 row has all zeros, delete this column."

drawcorplot(cor_no_cb2[,grep('Site3.', colnames(cor_no_cb2))])
## [1] "24 row has all zeros, delete this column."

drawcorplot(cor_no_cb2[,grep('Site4.', colnames(cor_no_cb2))])
## [1] "18 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."

drawcorplot(cor_no_cb2[,grep('Site5.', colnames(cor_no_cb2))])
## [1] "17 row has all zeros, delete this column."
## [1] "19 row has all zeros, delete this column."
## [1] "22 row has all zeros, delete this column."
## [1] "22 row has all zeros, delete this column."

combat

drawcorplot((cor_cb2[,grep('Site1.', colnames(cor_cb2))]))
## [1] "20 row has all zeros, delete this column."

drawcorplot((cor_cb2[,grep('Site2.', colnames(cor_cb2))]))
## [1] "17 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "30 row has all zeros, delete this column."

drawcorplot((cor_cb2[,grep('Site3.', colnames(cor_cb2))]))
## [1] "15 row has all zeros, delete this column."

drawcorplot((cor_cb2[,grep('Site4.', colnames(cor_cb2))]))
## [1] "8 row has all zeros, delete this column."
## [1] "21 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."

drawcorplot((cor_cb2[,grep('Site5.', colnames(cor_cb2))]))
## [1] "17 row has all zeros, delete this column."
## [1] "18 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."

Agreement computation

  • top 20 component agreement (spatial correlation matching)

  • 1- sqrt(average( (I - correlation matrix)^2))

computemat <- function(mat,ncomp=20){
  tmp=rerow(mat)[1:ncomp,1:ncomp]
  re=data.frame(diag.element = mean(diag(tmp)), offdiag.element = 1-mean(tmp[lower.tri(tmp, diag=FALSE)])  )
  return(re)
}

dat=do.call(rbind, lapply(c(5,10,15,20,25),
                          function(j){
                            data.frame(site = 1:5, ncomp=j,
           combat = do.call(rbind,lapply(paste('Site',1:5,'.',sep=''),function(str)computemat(t(cor_cb[,grep(str, colnames(cor_cb))]),ncomp=j))),
           no.com=do.call(rbind,lapply(paste('Site',1:5,'.',sep=''),function(str)computemat(t(cor_no_cb[,grep(str, colnames(cor_no_cb))]),ncomp=j))))
                          }))
fig1=dat %>% 
  melt(., id=c('site','ncomp')) %>%
  mutate(method = factor(substr(variable,1,6), levels=c('combat','no.com'), labels=c('ComBat','no-ComBat'))) %>%
  mutate(element = factor(substr(variable,8,11), levels=c('diag','offd'), labels=c('mean(Diagonal Elements)','1-mean(Off-Diagnonal Elements)'))) %>%
  ggplot(., aes(x=site, y=value, color=method, group=method)) + 
  geom_point() + facet_grid(element~ncomp) +
  geom_line() + ylim(c(0.5,1)) +
  ylab('Similarity') + 
  theme_bw() +
    theme(legend.position='bottom')+ 
  labs(caption='* Higher values indicates higher similarity between matched components.\n * The similarity measures range from 0 to 1.')

fig1

png('OCDglobal_n140_SiteSimilarity_B1000.png', width=2400, height=1800, res=300)
fig1
dev.off()
## quartz_off_screen 
##                 2
fig1_1=dat %>% 
  select(-combat.offdiag.element,-no.com.offdiag.element) %>%
  melt(., id=c('site','ncomp')) %>%
  mutate(method = factor(substr(variable,1,6), levels=c('combat','no.com'), labels=c('ComBat','no-ComBat'))) %>%
  ggplot(., aes(x=site, y=value, color=method, group=method)) + 
  geom_point() + facet_grid(~ncomp) +
  geom_line() + ylim(c(0.5,1)) +
  ylab('Similarity') + 
  theme_bw() +
    theme(legend.position='bottom')+ 
  labs(caption='* Higher values indicates higher similarity between matched components.\n * The similarity measures range from 0 to 1.')

fig1_1

png('OCDglobal_n140_SiteSimilarity_diag_B1000.png', width=2400, height=1200, res=300)
fig1_1
dev.off()
## quartz_off_screen 
##                 2

Combat vs. no combat comparison

combatmat=abs(cor(all_cb,all_nocb))
combatmat = combatmat[, apply(combatmat, 2, function(s)all(is.na(s))) == FALSE]
combatmat = combatmat[apply(combatmat, 1, function(s)all(is.na(s))) == FALSE, ]


thresh=0.5
combatmat2=combatmat;combatmat2[abs(combatmat)<thresh]<-0
combatmat2 = combatmat2[apply(combatmat2,1,function(ss)all(ss<thresh))==FALSE, apply(combatmat2,2,function(ss)all(ss< thresh))==FALSE]

#combatmat2 = combatmat2[apply(combatmat2,1,function(ss)all(ss<thresh))==FALSE, apply(combatmat2,2,function(ss)all(ss< thresh))==FALSE]

#combatmat2 = combatmat2[, apply(combatmat2,2,function(ss)all(ss<thresh))==FALSE]
#corrplot(rerow(combatmat2),col = colist, col.lim=c(0,1),is.corr = FALSE)
 
#corrplot(rerow(t(rerow(combatmat2))),col = COL1('YlOrBr', 200), col.lim=c(0,1),is.corr = FALSE)
#corrplot((rerow(combatmat2)),col = COL1('YlOrBr', 200), col.lim=c(0,1),is.corr = FALSE)

thresh=0.5
rerowed.mat = rerow(t(combatmat2))[c(1:31,33,32),]
## [1] "8 row has all zeros, delete this column."
corrplot(t(rerowed.mat),col = colist, col.lim=c(0,1),is.corr = FALSE)

pp=min(nrow(rerowed.mat),ncol(rerowed.mat));
matched.dat = data.frame(combatics = colnames(rerowed.mat)[1:pp],nocombatics = row.names(rerowed.mat)[1:pp], 
                         matchingcor = diag(rerowed.mat[1:pp,1:pp])) %>%
  filter(matchingcor>0.5)
knitr::kable(matched.dat) %>% kableExtra::kable_classic(full=F)
combatics nocombatics matchingcor
Combat,All.IC01 noCombat.All.IC01 0.9954084
Combat,All.IC02 noCombat.All.IC02 0.9987395
Combat,All.IC03 noCombat.All.IC03 0.9875149
Combat,All.IC04 noCombat.All.IC04 0.9967404
Combat,All.IC05 noCombat.All.IC06 0.9919961
Combat,All.IC06 noCombat.All.IC05 0.8671488
Combat,All.IC07 noCombat.All.IC09 0.9887035
Combat,All.IC09 noCombat.All.IC17 0.7780303
Combat,All.IC10 noCombat.All.IC12 0.9805605
Combat,All.IC11 noCombat.All.IC11 0.9772330
Combat,All.IC12 noCombat.All.IC25 0.9079190
Combat,All.IC13 noCombat.All.IC10 0.9648221
Combat,All.IC14 noCombat.All.IC13 0.9941005
Combat,All.IC15 noCombat.All.IC18 0.9505280
Combat,All.IC16 noCombat.All.IC20 0.9584342
Combat,All.IC17 noCombat.All.IC16 0.9743688
Combat,All.IC18 noCombat.All.IC22 0.9728010
Combat,All.IC19 noCombat.All.IC26 0.9731191
Combat,All.IC20 noCombat.All.IC30 0.6589200
Combat,All.IC21 noCombat.All.IC23 0.9628430
Combat,All.IC22 noCombat.All.IC24 0.9302998
Combat,All.IC23 noCombat.All.IC15 0.8817485
Combat,All.IC24 noCombat.All.IC14 0.9196102
Combat,All.IC25 noCombat.All.IC21 0.9062225
Combat,All.IC26 noCombat.All.IC29 0.8700443
Combat,All.IC27 noCombat.All.IC19 0.8775439
Combat,All.IC28 noCombat.All.IC28 0.9717423
Combat,All.IC29 noCombat.All.IC32 0.8331042
Combat,All.IC30 noCombat.All.IC33 0.9382460
Combat,All.IC31 noCombat.All.IC31 0.9338284
Combat,All.IC32 noCombat.All.IC35 0.9173853
Combat,All.IC33 noCombat.All.IC34 0.9226888
Combat,All.IC34 noCombat.All.IC08 0.5759699 
#corrplot((rerow((combatmat2))),col = colist, col.lim=c(0,1),is.corr = FALSE)

#corrplot(rerow(t(rerow((combatmat2)))),col = colist, col.lim=c(0,1),is.corr = FALSE)

#corrplot(rerow(combatmat2),col = colist, col.lim=c(0,1),is.corr = FALSE)

Explore site stability of LICs with ComBat

matchsite <- function(mat, threshold=0.5){
  mat[abs(mat)<threshold]<-0
  varnames=rownames(mat)
  mat = mat[apply(mat, 1, function(s)all(abs(s<threshold)))==FALSE,]
  mat = mat[,apply(mat, 2, function(s)all(abs(s<threshold)))==FALSE]
  rerowed.mat = rerow(t(mat))
  pp=min(nrow(rerowed.mat),ncol(rerowed.mat));
  tmpdat=data.frame(col1=varnames) %>%
    left_join(.,data.frame(col1=colnames(rerowed.mat)[1:pp],col2=row.names(rerowed.mat)[1:pp], col3=diag(rerowed.mat[1:pp,1:pp])))
  tmpdat$col3[ tmpdat$col3==0]<-NA
  return(tmpdat)
}

figs = lapply(1:5,function(jj)matchsite((cor_cb[,grep(paste('Site',jj,'.',sep=''), colnames(cor_cb))])))
## [1] "20 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "30 row has all zeros, delete this column."
## [1] "15 row has all zeros, delete this column."
## [1] "8 row has all zeros, delete this column."
## [1] "21 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "18 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
combat.sitematch = figs[[1]] %>% rename(site1=col2, site1.r=col3) %>% 
  left_join(.,figs[[2]] %>% rename(site2=col2, site2.r=col3)) %>% 
  left_join(.,figs[[3]] %>% rename(site3=col2, site3.r=col3))%>% 
  left_join(.,figs[[4]] %>% rename(site4=col2, site4.r=col3))%>% 
  left_join(.,figs[[5]] %>% rename(site5=col2, site5.r=col3)) %>%
  rename(combatics = col1)


combat.sitematch_cor = combat.sitematch %>%
  dplyr::select(combatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
  left_join(., 
            combat.sitematch %>%
              dplyr::select(combatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
              reshape2::melt(., id='combatics') %>% group_by(combatics) %>%
              summarise(num.matched = 5-sum(is.na(value)), mean.r = mean(value, na.rm=TRUE))
  )

combat.sitematch_cor %>%
  knitr::kable(., digits=3,  caption='ComBat, Site stability') %>%
  kableExtra::kable_classic(.,full=F)
ComBat, Site stability
combatics site1.r site2.r site3.r site4.r site5.r num.matched mean.r
Combat,All.IC01 0.998 0.998 0.996 0.998 0.997 5 0.997
Combat,All.IC02 0.999 0.999 0.999 0.999 0.998 5 0.999
Combat,All.IC03 0.990 0.989 0.994 0.984 0.989 5 0.989
Combat,All.IC04 0.973 0.907 0.944 0.989 0.981 5 0.959
Combat,All.IC05 0.989 0.990 0.993 0.984 0.984 5 0.988
Combat,All.IC06 0.937 0.771 0.919 0.974 0.931 5 0.907
Combat,All.IC07 0.975 0.951 0.978 0.929 0.982 5 0.963
Combat,All.IC08 0.560 0.962 0.938 NA 0.543 4 0.751
Combat,All.IC09 0.899 0.572 0.570 0.576 0.813 5 0.686
Combat,All.IC10 0.923 0.902 0.885 0.903 0.875 5 0.898
Combat,All.IC11 0.852 0.930 0.924 0.972 0.709 5 0.877
Combat,All.IC12 0.579 0.926 0.931 0.812 0.730 5 0.796
Combat,All.IC13 0.996 0.995 0.933 0.961 0.894 5 0.956
Combat,All.IC14 0.943 0.800 0.758 0.935 0.897 5 0.867
Combat,All.IC15 0.907 0.780 NA 0.839 0.841 4 0.842
Combat,All.IC16 0.810 0.661 0.972 0.783 0.723 5 0.790
Combat,All.IC17 0.595 NA NA NA NA 1 0.595
Combat,All.IC18 NA 0.914 0.867 0.967 0.961 4 0.927
Combat,All.IC19 0.778 0.577 0.828 0.922 NA 4 0.776
Combat,All.IC20 0.906 0.528 0.705 NA 0.869 4 0.752
Combat,All.IC21 NA NA 0.829 0.842 0.587 3 0.753
Combat,All.IC22 0.912 0.927 0.926 NA 0.958 4 0.931
Combat,All.IC23 0.901 0.742 0.616 0.752 0.936 5 0.789
Combat,All.IC24 0.960 0.813 0.948 0.873 0.829 5 0.885
Combat,All.IC25 NA 0.735 0.777 NA NA 2 0.756
Combat,All.IC26 NA NA NA 0.775 NA 1 0.775
Combat,All.IC27 0.915 0.892 0.993 0.689 0.986 5 0.895
Combat,All.IC28 0.962 0.930 0.934 0.514 0.590 5 0.786
Combat,All.IC29 0.861 0.931 0.814 NA 0.585 4 0.798
Combat,All.IC30 0.582 0.661 NA 0.883 0.512 4 0.659
Combat,All.IC31 NA NA NA NA 0.888 1 0.888
Combat,All.IC32 0.897 0.939 0.930 0.899 0.894 5 0.912
Combat,All.IC33 0.913 0.957 0.844 0.926 0.919 5 0.912
Combat,All.IC34 0.754 NA 0.942 0.746 0.931 4 0.843
Combat,All.IC35 0.766 0.507 0.648 NA NA 3 0.640
Combat,All.IC36 NA NA 0.712 NA NA 1 0.712
Combat,All.IC37 NA NA 0.917 NA NA 1 0.917

Explore site stability of LICs without ComBat

figs_nocombat = lapply(1:5,function(jj)matchsite((cor_no_cb[,grep(paste('Site',jj,'.',sep=''), colnames(cor_no_cb))])))
## [1] "16 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "18 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "19 row has all zeros, delete this column."
## [1] "22 row has all zeros, delete this column."
## [1] "22 row has all zeros, delete this column."
nocombat.sitematch = figs_nocombat[[1]] %>% rename(site1=col2, site1.r=col3) %>% 
  left_join(.,figs_nocombat[[2]] %>% rename(site2=col2, site2.r=col3)) %>% 
  left_join(.,figs_nocombat[[3]] %>% rename(site3=col2, site3.r=col3))%>% 
  left_join(.,figs_nocombat[[4]] %>% rename(site4=col2, site4.r=col3))%>% 
  left_join(.,figs_nocombat[[5]] %>% rename(site5=col2, site5.r=col3)) %>%
  rename(nocombatics = col1)

npcombat.sitematch_cor = nocombat.sitematch %>%
  dplyr::select(nocombatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
  left_join(., 
            nocombat.sitematch %>%
              dplyr::select(nocombatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
              reshape2::melt(., id='nocombatics') %>% group_by(nocombatics) %>%
              summarise(num.matched = 5-sum(is.na(value)), mean.r = mean(value, na.rm=TRUE))
  )

npcombat.sitematch_cor %>%
  knitr::kable(., digits=3,  caption='noComBat, Site stability') %>%
  kableExtra::kable_classic(.,full=F)
noComBat, Site stability
nocombatics site1.r site2.r site3.r site4.r site5.r num.matched mean.r
noCombat.All.IC01 0.998 0.997 0.992 0.995 0.994 5 0.995
noCombat.All.IC02 0.999 0.999 0.999 0.999 0.999 5 0.999
noCombat.All.IC03 0.996 0.986 0.991 0.994 0.993 5 0.992
noCombat.All.IC04 0.986 0.903 0.937 0.926 0.984 5 0.947
noCombat.All.IC05 0.843 0.953 0.812 0.934 0.922 5 0.893
noCombat.All.IC06 0.992 0.966 0.992 0.925 0.988 5 0.972
noCombat.All.IC07 0.944 0.940 0.798 0.804 0.763 5 0.850
noCombat.All.IC08 0.875 NA 0.870 0.638 0.913 4 0.824
noCombat.All.IC09 0.977 0.903 0.925 0.972 0.977 5 0.951
noCombat.All.IC10 0.980 0.970 0.951 0.952 0.893 5 0.949
noCombat.All.IC11 0.901 0.934 0.912 0.893 0.726 5 0.873
noCombat.All.IC12 0.802 0.855 0.583 0.869 0.868 5 0.795
noCombat.All.IC13 0.943 0.897 0.881 0.947 0.898 5 0.913
noCombat.All.IC14 0.961 0.938 0.869 0.834 NA 4 0.900
noCombat.All.IC15 0.977 NA 0.844 0.964 0.940 4 0.932
noCombat.All.IC16 NA NA 0.747 0.554 0.509 3 0.603
noCombat.All.IC17 NA 0.539 0.543 0.634 0.611 4 0.582
noCombat.All.IC18 0.892 0.930 0.838 NA NA 3 0.887
noCombat.All.IC19 0.991 0.757 0.951 0.872 0.924 5 0.899
noCombat.All.IC20 0.958 0.628 0.961 0.539 0.609 5 0.739
noCombat.All.IC21 0.707 0.908 0.746 NA NA 3 0.787
noCombat.All.IC22 0.666 0.719 0.779 0.940 0.969 5 0.815
noCombat.All.IC23 0.626 0.754 0.835 0.843 0.775 5 0.767
noCombat.All.IC24 0.834 0.821 NA NA 0.877 3 0.844
noCombat.All.IC25 NA 0.701 NA 0.916 NA 2 0.808
noCombat.All.IC26 NA NA 0.800 0.907 NA 2 0.853
noCombat.All.IC27 NA NA 0.634 NA 0.898 2 0.766
noCombat.All.IC28 0.965 NA NA NA NA 1 0.965
noCombat.All.IC29 NA NA NA NA NA 0 NaN
noCombat.All.IC30 0.565 NA NA NA NA 1 0.565
noCombat.All.IC31 NA NA NA NA NA 0 NaN
noCombat.All.IC32 0.961 0.821 0.945 NA 0.789 4 0.879
noCombat.All.IC33 0.641 NA NA 0.788 NA 2 0.715
noCombat.All.IC34 NA 0.878 NA NA NA 1 0.878
noCombat.All.IC35 NA NA 0.862 NA NA 1 0.862 
aa = nocombat.sitematch %>% 
  dplyr::select(nocombatics, site1, site2, site3, site4, site5) %>%
  left_join(., npcombat.sitematch_cor %>% select(nocombatics, num.matched,mean.r)) %>%
    dplyr::filter(num.matched == 5 & mean.r>0.8) 

figs = lapply(1:5,
              function(j){
                obj0 = site$loso[[j]][[1]]
                obj=obj0[[3]]
   
                tmp=do.call(cbind,lapply(obj, unlist)) %>% 
                  data.frame(.) %>% 
                  set_names(c('CT','AREA','FAb1000','MDb1000','FC')) 
                W40 = tmp[order(site$loso[[j]][[1]][[15]], decreasing=TRUE),] %>%
                  mutate(IC=paste('Site',j,'.IC',stringr::str_pad(1:nrow(tmp),2, pad = "0"), sep=''))
                W40[,1:5]=abs(W40[,1:5])/apply(abs(W40[,1:5]), 1,sum)
                row.names(W40)<-W40$IC
                
                fig = W40[aa[,j+1],] %>%
                  mutate(IC=1:nrow(aa)) %>% 
                  melt(., id='IC') %>%
                  ggplot(.,aes(x=reorder(IC, desc(IC)),y=value,fill=variable))+
                  geom_bar(stat = "identity",position="fill", color="black", width=0.9) + 
                  ggtitle(paste('Site',j)) +
                  ylab("Modality Contribution (%)") + coord_flip() + scale_y_continuous(labels=scales::percent) +
                  theme_minimal() + scale_fill_manual(values=c('salmon','salmon3','aquamarine1',
                                               'deepskyblue1','gold')) + 
                  theme(legend.position='none',
                        axis.title.y=element_blank(),
                        axis.text.y=element_blank(),
                        axis.ticks.y=element_blank())
                return(list(fig,W40[aa[,j+1],]))
              })

Morig40 = all$Morig40
W40=data.frame(do.call(cbind,lapply(Morig40[[3]], unlist)))[order(all$Morig40[[15]], decreasing=TRUE),] %>% 
  set_names(c('CT','AREA','FAB1K','MDB1K','FC')) 
W40=W40%>%
  mutate(IC=paste('noCombat.All.IC',stringr::str_pad(1:nrow(W40),2, pad = "0"), sep=''))
W40[,1:5]=abs(W40[,1:5])/apply(abs(W40[,1:5]), 1,sum)
row.names(W40)<-W40$IC
fig.all = W40[aa[,1],] %>%
                  melt(., id='IC') %>%
  mutate(IC = gsub('noCombat.All.','',IC)) %>%
                  ggplot(.,aes(x=reorder(IC, desc(IC)),y=value,fill=variable))+
                  geom_bar(stat = "identity",position="fill", color="black", width=0.9) + 
                  ggtitle(paste('All Sites')) +
                  xlab('LICs before ComBat')+
                  ylab("Modality Contribution (%)") + coord_flip() + scale_y_continuous(labels=scales::percent) +
                  theme_minimal() + scale_fill_manual(values=c('salmon','salmon3','aquamarine1',
                                               'deepskyblue1','gold')) + 
                  theme(legend.position='none')

grid.arrange(fig.all, figs[[1]][[1]], figs[[2]][[1]], figs[[3]][[1]], figs[[4]][[1]], figs[[5]][[1]], 
             widths=c(1.5,1,1,1,1,1))

#knitr::kable(rbind(figs[[5]][[2]],W40[aa[,1],]),digits=2)
dims = do.call(rbind,lapply(Morig40[[7]],function(obj)dim(obj[[1]])))[,1]

Explore site stability of LICs after ComBat

figs_combat = lapply(1:5,function(jj)matchsite((cor_cb[,grep(paste('Site',jj,'.',sep=''), colnames(cor_cb))])))
## [1] "20 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "30 row has all zeros, delete this column."
## [1] "15 row has all zeros, delete this column."
## [1] "8 row has all zeros, delete this column."
## [1] "21 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "18 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
combat.sitematch = figs_combat[[1]] %>% rename(site1=col2, site1.r=col3) %>% 
  left_join(.,figs_combat[[2]] %>% rename(site2=col2, site2.r=col3)) %>% 
  left_join(.,figs_combat[[3]] %>% rename(site3=col2, site3.r=col3))%>% 
  left_join(.,figs_combat[[4]] %>% rename(site4=col2, site4.r=col3))%>% 
  left_join(.,figs_combat[[5]] %>% rename(site5=col2, site5.r=col3)) %>%
  rename(combatics = col1)

combat.sitematch_cor = combat.sitematch %>%
  dplyr::select(combatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
  left_join(., 
            combat.sitematch %>%
              dplyr::select(combatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
              reshape2::melt(., id='combatics') %>% group_by(combatics) %>%
              summarise(mean.r = mean(value, na.rm=TRUE),num.matched = 5-sum(is.na(value)))
  )

combat.sitematch_cor %>%
  knitr::kable(., digits=3,  caption='ComBat, Site stability') %>%
  kableExtra::kable_classic(.,full=F)
ComBat, Site stability
combatics site1.r site2.r site3.r site4.r site5.r mean.r num.matched
Combat,All.IC01 0.998 0.998 0.996 0.998 0.997 0.997 5
Combat,All.IC02 0.999 0.999 0.999 0.999 0.998 0.999 5
Combat,All.IC03 0.990 0.989 0.994 0.984 0.989 0.989 5
Combat,All.IC04 0.973 0.907 0.944 0.989 0.981 0.959 5
Combat,All.IC05 0.989 0.990 0.993 0.984 0.984 0.988 5
Combat,All.IC06 0.937 0.771 0.919 0.974 0.931 0.907 5
Combat,All.IC07 0.975 0.951 0.978 0.929 0.982 0.963 5
Combat,All.IC08 0.560 0.962 0.938 NA 0.543 0.751 4
Combat,All.IC09 0.899 0.572 0.570 0.576 0.813 0.686 5
Combat,All.IC10 0.923 0.902 0.885 0.903 0.875 0.898 5
Combat,All.IC11 0.852 0.930 0.924 0.972 0.709 0.877 5
Combat,All.IC12 0.579 0.926 0.931 0.812 0.730 0.796 5
Combat,All.IC13 0.996 0.995 0.933 0.961 0.894 0.956 5
Combat,All.IC14 0.943 0.800 0.758 0.935 0.897 0.867 5
Combat,All.IC15 0.907 0.780 NA 0.839 0.841 0.842 4
Combat,All.IC16 0.810 0.661 0.972 0.783 0.723 0.790 5
Combat,All.IC17 0.595 NA NA NA NA 0.595 1
Combat,All.IC18 NA 0.914 0.867 0.967 0.961 0.927 4
Combat,All.IC19 0.778 0.577 0.828 0.922 NA 0.776 4
Combat,All.IC20 0.906 0.528 0.705 NA 0.869 0.752 4
Combat,All.IC21 NA NA 0.829 0.842 0.587 0.753 3
Combat,All.IC22 0.912 0.927 0.926 NA 0.958 0.931 4
Combat,All.IC23 0.901 0.742 0.616 0.752 0.936 0.789 5
Combat,All.IC24 0.960 0.813 0.948 0.873 0.829 0.885 5
Combat,All.IC25 NA 0.735 0.777 NA NA 0.756 2
Combat,All.IC26 NA NA NA 0.775 NA 0.775 1
Combat,All.IC27 0.915 0.892 0.993 0.689 0.986 0.895 5
Combat,All.IC28 0.962 0.930 0.934 0.514 0.590 0.786 5
Combat,All.IC29 0.861 0.931 0.814 NA 0.585 0.798 4
Combat,All.IC30 0.582 0.661 NA 0.883 0.512 0.659 4
Combat,All.IC31 NA NA NA NA 0.888 0.888 1
Combat,All.IC32 0.897 0.939 0.930 0.899 0.894 0.912 5
Combat,All.IC33 0.913 0.957 0.844 0.926 0.919 0.912 5
Combat,All.IC34 0.754 NA 0.942 0.746 0.931 0.843 4
Combat,All.IC35 0.766 0.507 0.648 NA NA 0.640 3
Combat,All.IC36 NA NA 0.712 NA NA 0.712 1
Combat,All.IC37 NA NA 0.917 NA NA 0.917 1 
aa = combat.sitematch %>% 
  dplyr::select(combatics, site1, site2, site3, site4, site5) %>%
  left_join(., combat.sitematch_cor %>% select(combatics, num.matched,mean.r)) %>%
    dplyr::filter(num.matched==5 & mean.r>0.8)
                  
figs = lapply(1:5,
              function(j){
                obj0 = site$loso.comba[[j]][[1]]
                obj=obj0[[3]]
   
                tmp=do.call(cbind,lapply(obj, unlist)) %>% 
                  data.frame(.) %>% 
                  set_names(c('CT','AREA','FAb1000','MDb1000','FC')) 
                W40 = tmp[order(site$loso.comba[[j]][[1]][[15]], decreasing=TRUE),] %>%
                  mutate(IC=paste('Site',j,'.IC',stringr::str_pad(1:nrow(tmp),2, pad = "0"), sep=''))
                W40[,1:5]=abs(W40[,1:5])/apply(abs(W40[,1:5]), 1,sum)
                row.names(W40)<-W40$IC
                
                fig = W40[aa[,j+1],] %>%
                  mutate(IC=1:nrow(aa)) %>% 
                  melt(., id='IC') %>%
                  ggplot(.,aes(x=reorder(IC, desc(IC)),y=value,fill=variable))+
                  geom_bar(stat = "identity",position="fill", color="black", width=0.9) + 
                  ggtitle(paste('Site',j)) +
                  ylab("Modality Contribution (%)") + coord_flip() + scale_y_continuous(labels=scales::percent) +
                  theme_minimal() + scale_fill_manual(values=c('salmon','salmon3','aquamarine1',
                                               'deepskyblue1','gold')) + 
                  theme(legend.position='none',
                        axis.title.y=element_blank(),
                        axis.text.y=element_blank(),
                        axis.ticks.y=element_blank())
                return(list(fig,W40[aa[,j+1],]))
              })

Morig40 = all$Morig40.combat
W40=data.frame(do.call(cbind,lapply(Morig40[[3]], unlist)))[order(Morig40[[15]], decreasing=TRUE),] %>% 
  set_names(c('CT','AREA','FAB1K','MDB1K','FC')) 
W40=W40%>%
  mutate(IC=paste('Combat,All.IC',stringr::str_pad(1:nrow(W40),2, pad = "0"), sep=''))
W40[,1:5]=abs(W40[,1:5])/apply(abs(W40[,1:5]), 1,sum)
row.names(W40)<-W40$IC
fig.all = W40[aa[,1],] %>%
                  melt(., id='IC') %>%
  mutate(IC = gsub('Combat.All.','',IC)) %>%
                  ggplot(.,aes(x=reorder(IC, desc(IC)),y=value,fill=variable))+
                  geom_bar(stat = "identity",position="fill", color="black", width=0.9) + 
                  ggtitle(paste('All Sites')) +
                  xlab('LICs after ComBat')+
                  ylab("Modality Contribution (%)") + coord_flip() + scale_y_continuous(labels=scales::percent) +
                  theme_minimal() + scale_fill_manual(values=c('salmon','salmon3','aquamarine1',
                                               'deepskyblue1','gold')) + 
                  theme(legend.position='none')

grid.arrange(fig.all, figs[[1]][[1]], figs[[2]][[1]], figs[[3]][[1]], figs[[4]][[1]], figs[[5]][[1]], 
             widths=c(1.5,1,1,1,1,1))